Filter criteria:
Filter differentially expressed genes between autism and control (adj. p-value < 0.05 and lfc<log2(1.2))
No samples are removed based on network connectivity z-scores
# # FILTER DE GENES
# load('./../working_data/RNAseq_ASD_4region_normalized_vst.Rdata')
#
# # Balance Groups by covariates, remove singular batches (none)
# to_keep = (datMeta$Subject_ID != 'AN03345') & !is.na(datMeta$Dx)
# table(to_keep)
# datMeta = datMeta[to_keep,]
# datExpr = datExpr[,to_keep]
#
# # Select genes differentially expressed in ASD
# mod = model.matrix(~ Dx, data=datMeta)
# corfit = duplicateCorrelation(datExpr, mod, block=datMeta$Subject_ID)
# lmfit = lmFit(datExpr, design=mod, block=datMeta$Subject_ID, correlation=corfit$consensus)
#
# fit = eBayes(lmfit, trend=T, robust=T)
# top_genes = topTable(fit, coef=2, number=nrow(datExpr))
# ordered_top_genes = top_genes[match(rownames(datExpr), rownames(top_genes)),]
#
# ASD_pvals = rownames(ordered_top_genes)[ordered_top_genes$adj.P.Val<0.05] # 6% of genes (3380)
# ASD_lfc = rownames(ordered_top_genes)[abs(ordered_top_genes$logFC)>log2(1.2)] # 13% of genes (6732)
# ASD_pvals_lfc = intersect(ASD_pvals, ASD_lfc) # 6% of genes (2928)
#
# datExpr = datExpr[match(ASD_pvals_lfc, rownames(datExpr)),]
#
# datProbes = datProbes[rownames(datProbes) %in% rownames(datExpr),]
#
# rm(to_keep, corfit, fit, lmfit, mod, ASD_pvals, ASD_lfc, ASD_pvals_lfc, top_genes)
#
# save(file='./../working_data/RNAseq_ASD_4region_DEgenes_vst_adj_pval_lfc.Rdata', datExpr, datMeta, datProbes, ordered_top_genes)
load('./../working_data/RNAseq_ASD_4region_DEgenes_vst_adj_pval_lfc.Rdata')
DE_info = ordered_top_genes
rm(ordered_top_genes)
glue('Number of genes: ', nrow(datExpr), '\n',
'Number of samples: ', ncol(datExpr), ' (', sum(datMeta$Diagnosis_=='ASD'), ' ASD, ',
sum(datMeta$Diagnosis_!='ASD'), ' controls)')
## Number of genes: 3071
## Number of samples: 86 (51 ASD, 35 controls)
reduce_dim_datExpr = function(datExpr, datMeta, var_explained=0.8, filter_controls=FALSE){
datExpr = data.frame(datExpr)
if(filter_controls){
datMeta = datMeta %>% filter(Diagnosis_=='ASD')
datExpr = datExpr %>% select(paste0('X', datMeta_ASD$Dissected_Sample_ID))
}
datExpr_pca = prcomp(t(datExpr), scale=TRUE)
last_pc = data.frame(summary(datExpr_pca)$importance[3,]) %>% rownames_to_column(var='ID') %>%
filter(.[[2]] >= var_explained) %>% top_n(-1) %>% dplyr::select(ID)
par(mfrow=c(1,2))
plot(summary(datExpr_pca)$importance[2,], type='b')
abline(v=substr(last_pc$ID, 3, nchar(last_pc$ID)), col='blue')
plot(summary(datExpr_pca)$importance[3,], type='b')
abline(h=var_explained, col='blue')
print(glue('Keeping top ', substr(last_pc$ID, 3, nchar(last_pc$ID)), ' components that explain ',
var_explained*100, '% of the variance'))
datExpr_top_pc = datExpr_pca$x %>% data.frame %>% dplyr::select(PC1:last_pc$ID)
return(list('datExpr'=datExpr_top_pc, 'datMeta'=datMeta, 'pca_output'=datExpr_pca))
}
reduce_dim_output = reduce_dim_datExpr(datExpr, datMeta)
## Keeping top 10 components that explain 80% of the variance
datExpr_redDim = reduce_dim_output$datExpr
datMeta_redDim = reduce_dim_output$datMeta
pca_output = reduce_dim_output$pca_output
rm(datSeq, datProbes, reduce_dim_datExpr, reduce_dim_output, datExpr, datMeta)
clusterings = list()
No recognisable best k, so chose k=5
set.seed(123)
wss = sapply(1:15, function(k) kmeans(datExpr_redDim, k, nstart=25)$tot.withinss)
plot(wss, type='b', main='K-Means Clustering')
best_k = 3
abline(v = best_k, col='blue')
datExpr_k_means = kmeans(datExpr_redDim, best_k, nstart=25)
clusterings[['km']] = datExpr_k_means$cluster
Chose k=7 as best number of clusters.
Clusters seem to be able to separate ASD and control samples pretty well and there are no noticeable patterns regarding sex, age or brain region in any cluster.
Younger ASD samples seem to be more similar to control samples than older ASD samples (pink cluster has most of the youngest samples), perhaps oder samples were only diagnosed in extreme cases and now milder ASD cases get diagnosed a well and milder cases have milder gene expression differences? The yellow cluster is made of young ASD samples. Most old ASD samples are also close together.
Colors:
Diagnosis: Blue=control, Green=ASD
Sex: Pink=Female, Blue=Male
Brain region: Pink=Frontal, Green=Temporal, Blue=Parietal, Purple=Occipital
Age: Purple=youngest, Yellow=oldest
h_clusts = datExpr_redDim %>% dist %>% hclust %>% as.dendrogram
# h_clusts %>% plot
best_k = 7
clusterings[['hc']] = cutree(h_clusts, best_k)
create_viridis_dict = function(){
min_age = datMeta_redDim$Age %>% min
max_age = datMeta_redDim$Age %>% max
viridis_age_cols = viridis(max_age - min_age + 1)
names(viridis_age_cols) = seq(min_age, max_age)
return(viridis_age_cols)
}
viridis_age_cols = create_viridis_dict()
dend_meta = datMeta_redDim[match(gsub('X','',labels(h_clusts)), rownames(datMeta_redDim)),] %>%
mutate('Diagnosis' = ifelse(Diagnosis_=='CTL','#008080','#86b300'), # Blue control, Green ASD
'Sex' = ifelse(Sex=='F','#ff6666','#008ae6'), # Pink Female, Blue Male
'Region' = case_when(Brain_lobe=='Frontal'~'#F8766D', # ggplot defaults for 4 colours
Brain_lobe=='Temporal'~'#7CAE00',
Brain_lobe=='Parietal'~'#00BFC4',
Brain_lobe=='Occipital'~'#C77CFF'),
'Age' = viridis_age_cols[as.character(Age)]) %>% # Purple: young, Yellow: old
dplyr::select(Age, Region, Sex, Diagnosis)
h_clusts %>% set('labels', rep('', nrow(datMeta_redDim))) %>% set('branches_k_color', k=best_k) %>% plot
colored_bars(colors=dend_meta)
Samples are grouped into two big clusters with 8 and 7 subclusters, respectively.
*Output plots in clustering_samples_04_05 folder
Following this paper’s guidelines:
Run PCA and keep enough components to explain 60% of the variance
Run ICA with that same number of nbComp as principal components kept to then filter them
Select components with kurtosis > 3
Assign obs to genes with FDR<0.01 using the fdrtool package
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
##
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
Not a good method for clustering samples because:
ICA does not perform well with small samples (see Figure 4 of this paper)
Warnings: (Warning in fdrtool(x, plot = F): There may be too few input test statistics for reliable FDR calculations!)
Leaves almost half of the observations (59) without a cluster:
ICA_clusters %>% rowSums %>% table
## .
## 0 1 2 3
## 59 24 2 1
ICA_clusters %>% mutate(cl_sum=rowSums(.)) %>% as.matrix %>% melt %>% ggplot(aes(Var2,Var1)) +
geom_tile(aes(fill=value)) + xlab('Clusters') + ylab('Samples') +
theme_minimal() + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + coord_flip()
SFT.R.sq starts in 0.83 and then 0.846 but it then decreases and never achieves a higher value, so chose power=2
best_power = datExpr_redDim %>% t %>% pickSoftThreshold(powerVector = seq(1, 30, by=2))
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.8380 1.6600 0.812 43.400 45.600 58.30
## 2 3 0.0135 0.0343 -0.254 19.500 20.100 35.80
## 3 5 0.2460 -0.3180 0.149 11.300 10.700 25.30
## 4 7 0.6310 -0.5260 0.624 7.320 6.310 19.10
## 5 9 0.6590 -0.6920 0.616 5.110 3.930 15.10
## 6 11 0.6950 -0.7960 0.735 3.760 2.640 12.30
## 7 13 0.7230 -0.7810 0.761 2.870 1.770 10.20
## 8 15 0.8140 -0.8490 0.886 2.260 1.260 8.59
## 9 17 0.8430 -0.8920 0.885 1.820 0.962 7.35
## 10 19 0.7800 -0.9520 0.782 1.490 0.690 6.36
## 11 21 0.8040 -0.9870 0.799 1.240 0.510 5.55
## 12 23 0.7840 -1.0100 0.772 1.050 0.404 4.88
## 13 25 0.7840 -0.9680 0.774 0.898 0.333 4.31
## 14 27 0.8300 -0.9280 0.871 0.775 0.278 3.84
## 15 29 0.8090 -0.9530 0.811 0.676 0.234 3.43
network = datExpr_redDim %>% t %>% blockwiseModules(power=2, numericLabels=TRUE)
clusterings[['WGCNA']] = network$colors
names(clusterings[['WGCNA']]) = rownames(datExpr_redDim)
It finds two main clusters and leaves only 3 observations without a cluster
table(clusterings[['WGCNA']])
##
## 0 1 2
## 3 50 33
Points don’t seem to follow a Gaussian distribution no matter the number of clusters, chose 4 points following the best k from K-means because the methods are similar
n_clust = datExpr_redDim %>% Optimal_Clusters_GMM(max_clusters=80, criterion='BIC', plot_data=FALSE)
plot(n_clust, type='l', main='Bayesian Information Criterion to choose number of clusters')
best_k = 3 # copying k-means best_k
gmm = datExpr_redDim %>% GMM(best_k)
clusterings[['GMM']] = gmm$Log_likelihood %>% apply(1, which.max)
Plot of clusters with their centroids in gray
gmm_points = rbind(datExpr_redDim, setNames(data.frame(gmm$centroids), names(datExpr_redDim)))
gmm_labels = c(clusterings[['GMM']], rep(NA, best_k)) %>% as.factor
ggplotly(gmm_points %>% ggplot(aes(x=PC1, y=PC2, color=gmm_labels)) + geom_point() + theme_minimal())
rm(wss, datExpr_k_means, h_clusts, cc_output, cc_output_c1, cc_output_c2, best_k, ICA_output,
ICA_clusters_names, signals_w_kurtosis, n_clust, gmm, gmm_points, gmm_labels, network, dend_meta,
best_power, c, viridis_age_cols, create_viridis_dict)
Using Adjusted Rand Index:
All clusterings are pretty similar except for WGCNA and ICA
K-means and Hierarchical clustering are the most similar
ICA seems to cluster samples by age or by subject
Consensus clustering and GMM are the ones with the strongest relation to ASD, but the relation is not that strong
clusters_plus_phenotype = clusterings
clusters_plus_phenotype[['ICA_NA']] = is.na(clusters_plus_phenotype[['ICA_min']])
clusters_plus_phenotype[['Subject']] = datMeta_redDim$Subject_ID
clusters_plus_phenotype[['ASD']] = datMeta_redDim$Diagnosis_
clusters_plus_phenotype[['Region']] = datMeta_redDim$Brain_lobe
clusters_plus_phenotype[['Sex']] = datMeta_redDim$Sex
clusters_plus_phenotype[['Age']] = datMeta_redDim$Age
cluster_sim = data.frame(matrix(nrow = length(clusters_plus_phenotype), ncol = length(clusters_plus_phenotype)))
for(i in 1:(length(clusters_plus_phenotype))){
cluster1 = as.factor(clusters_plus_phenotype[[i]])
for(j in (i):length(clusters_plus_phenotype)){
cluster2 = as.factor(clusters_plus_phenotype[[j]])
cluster_sim[i,j] = adj.rand.index(cluster1, cluster2)
}
}
colnames(cluster_sim) = names(clusters_plus_phenotype)
rownames(cluster_sim) = colnames(cluster_sim)
cluster_sim = cluster_sim %>% as.matrix %>% round(2)
heatmap.2(x = cluster_sim, Rowv = FALSE, Colv = FALSE, dendrogram = 'none',
cellnote = cluster_sim, notecol = 'black', trace = 'none', key = FALSE,
cexRow = 1, cexCol = 1, margins = c(7,7))
rm(i, j, cluster1, cluster2, clusters_plus_phenotype, cluster_sim)
plot_points = datExpr_redDim %>% data.frame() %>% dplyr::select(PC1:PC3) %>%
mutate(ID = rownames(.), subject_ID = datMeta_redDim$Subject_ID,
km_clust = as.factor(clusterings[['km']]), hc_clust = as.factor(clusterings[['hc']]),
cc_l1_clust = as.factor(clusterings[['cc_l1']]), cc_clust = as.factor(clusterings[['cc_l2']]),
ica_clust = as.factor(clusterings[['ICA_min']]), n_ica_clust = as.factor(rowSums(ICA_clusters)),
gmm_clust = as.factor(clusterings[['GMM']]), wgcna_clust = as.factor(clusterings[['WGCNA']]),
sex = as.factor(datMeta_redDim$Sex), region = as.factor(datMeta_redDim$Brain_lobe),
diagnosis = as.factor(datMeta_redDim$Diagnosis_), age = datMeta_redDim$Age)
selectable_scatter_plot(plot_points, plot_points)